Structure of multi-component/multi- Yukawa 
mixtures 

L. Blumf § and M. Ariasf 

fDepartment of Physics, P.O. Box 23343, University of Puerto Rico, Rio Piedras, PR 
00931-3343 and Department of Mathematics, Hill Center Rutgers University, 
Piscataway N.J. 08854 

Abstract. Recent small angle scattering experiments reveal new peaks in the 
structure function S(k) of colloidal systems ( S.H. Chen et al [TJ), in a region that 
was inaccessible with older instruments. It has been increasingly evident that a single 
(or double) Yukawa MSA-closures cannot account for these observations, and three 
or more terms are needed. On the other hand the MSA is not sufficiently accurate0 
more accurate theories such as the HNC have been tried. But while the MSA is 
asymptotically exact at high densities jSjit does not satisfy the low density asymptotics. 
This has been corrected in the soft-MSA by adding exponential type terms. The 
results compared to experiment and simulation for liquid sodium by Rahman and 
Paskin ( as shown in 0j) are remarkably good. We use here a general closure of the 
Ornstein Zernike equation which is not necessarily the MSA closure [H] • 

M 

%(r) = £4V*" r A K^=K^6^S^; r > * iS (1) 

n=l 

with the boundary condition for gij(r) = for r < aij. This general closure of 
the Ornstein Zernike equation will go well beyond the MSA since it has been tested 
by Monte Carlo simulation for tetrahedral water 0, toroidal ion channels |Hj and 
polyelectrolytes For this closure we get for the Laplace transform of the pair 
correlation function an explicitly symmetric result 

f 1 1 M z X (m) f (m) 

D T (s) Is 2 s ^ s + z m 

(2) 

This function is also easily transformed into S(k) by replacing s => ik. For low density 
situations (dilute colloids) D T (s) ~ 1 + 0{p) and the S{k) is a sum of M Lorentzians. 
For hard sphere PY mixtures we get the simple (compare |10l 

where D T (s) is a scalar function. For polydisperse electrolytes in the MSA a simpler 
expression is also obtained (compare \12\ ) . An explicit continued fraction solution of 
the 1 component multiyukawa case is also given. 
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1. Introduction 

There are many problems of practical and academic interest that can be formulated as 
closures of some kind of either scalar or matrix Ornstein-Zernike (OZ) equation. These 
closures can always be expressed by a sum of exponentials, which do form a complete 
basis set if we allow for complex numbers fHJ ■ ^ is of practical interest to be able to 
relate small angle scattering experiments as directly as posible to theoretical (molecular) 
parameters. For a number of systems the MSA ([E1III3IIE]) has been generally adequate 
(PUED- The GOCM which is a the single peak-single Yukawa closure-MSA description 
is adequate in the simple cases. However, recent high resolution experiments have 
shown that to account for the additional peaks seen in the SANS experiments more 
yukawas are required 0IH1HH]. The MSA and HNC closures of the OZ equation slightly 
underestimate the height of the interaction peak in the structure factor, whereas a close 
comparison of the cluster peak position is somewhat hindered by the limited resolution in 
Q of the simulation data, because the largest simulation box is never big enough. In most 
part of the practical range of potential parameters, the theoretical predictions coming 
from the two closures have comparable accuracy. Therefore, for fitting the neutron 
scattering intensity distributions, analytical structure factors are preferable since the 
resolution of the equations is more stable and fast. 

Our present solution makes the direct comparison between theory and experiment 
feasible since a simple expression for the Fourier transform S(Q) is proposed which is 
directly related to the closure of the OZ equation. In fact with the very mild assumption 
made is that the direct correlation function can be expanded as 

M M 

Cij (r) = J2 Klfe-^-^/r = £ K$e"" r /r\ r > a t3 (3) 

n=l n=l 

In the factorizable case (which is also the electrostatic charge case), we have 

K\? = K^5^5l n) ; 4 n) = K^dfhf- 5^ = d^e^l 2 (4) 

with 

9ij(?) = f or r - °V 
Here z n can be a complex number, and then eq.© is a complete function basis set. 
This closure works very well with a number of potentials, that include Lennard- Jones, 
Buckingham, and liquid metals [311201 • The effective interaction among charged colloids 
and globular proteins in solution is also well described by this potential which is the 
sum of a short range part attraction, usually accounting for van der Waals or entropic 
forces, plus a long range repulsion (of up to a few particle diameters), arising from 
screened electrostatic potential[21j. The two- Yukawa (2Y) fluid, consisting of particles 
interacting with this potential can be usefully applied to a wide variety of systems such 
as C60 fullerenes at high temperature [22|, globular proteins |2HJ E3|, and the Derjaguin- 
Landau- Verwey-Overbeek (DLVO)^Hl potential for liophobic colloids. 
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While the initial motivation was to study simple approximations like the 
Mean Spherical (MSA) or Generalized Mean Spherical Approximation (GMSA), the 
availability of closed form solutions for the general closure of the hard core OZ equation 
makes it possible to write down analytical solutions for any given approximation that 
can be formulated by writing the direct correlation function c(r) outside the hard core 
as 

M M 

c (r) = Y, K {n) e- Zn{r - u) /r = Y IC {n) e- z " r /r (5) 

n=l n=l 

In this equation is the interaction/closure constant used in the general solution first 
found by Blum and Hoye (which we will call BH78) [20], while JC^ is the definition used 
in the later general solution by Blum, Vericat and Herrera (BVH92 in what follows) [26J. 
In this work we will use the more common notation of BVH92. The case of factored 
interactions was discussed by Blum, [23 EH] an d by Ginoza EOl EU E2j • The general 
solution of this problem is given in terms of the scaling matrix T which will comply 
with the physical constraints of the systems (Blum et al. ESI El])- The solution of 
the resulting algebraic equations has a number of branches (Pastore |35j). The rigorous 
analysis of this question is very complex. We have used the following working hypothesis: 

(i) The singularities of the equations are poles, branch cuts and essential ( exponential) 
sigularities. 

(ii) The physical branch is the one for which the correct zero density result is obtained. 
This is equivalent to what is done in the rigorous treatment van der Waals theory, 
and includes metastable regions and spinodals. 

(iii) The analysis of the repulsive part of the potential is more delicate since it is 
clearly related to the exponential singularities, and the convexity of the limiting 
high density configurations [SEj- 

For only one component the matrix T can be assumed to be diagonal without loss 
of generality, and explicit expressions for the closure relations for any arbitrary number 
of Yukawa exponents M were obtained. The solution is then remarkably simple in 
the MSA since then explicit formulas for the thermodynamic properties are obtained. 
However there are undetermined integration constants in the entropy which have to be 
adjusted to get the correct limiting behavior [T7] . 

2. Summary of Previous Work 

We study the Ornstein-Zernike (OZ) equation 

^(12) = Cij (12) + Y I d3h ik (13)p k c kj (32) (6) 
k J 

where hij (12) is the molecular total correlation function and Cy(12) is the molecular 
direct correlation function, pi is the number density of the molecules i, and % = 1,2 is 
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the position r*j , r 12 = \ f\ — r 2 | and is the distance of closest approach of two particles 
The direct correlation function is 

M 

Cij (r) = £ K§>e-*«<r-°^/r, r > a l3 (7) 

n=l 

and the pair correlation function is 

hij(r) = 9ij{r) - 1 = -1, r < ay (8) 
We use the Baxter- Wertheim (BW) factorization of the OZ equation 

[I + pH(k)][I-pC(A;)]=I (9) 
where I is the identity matrix, and we have used the notation 

roc 

H(k) = 2 drcos(kr)J(r) (10) 
Jo 

C(k) = 2 / dr cos(»S(r) (11) 
Jo 

The matrices J and S have matrix elements 

Jij{r) = 2tt / dsshijis) (12) 



Sij(r) = 2n J dsscij(s) (13) 

[I - pC(k)} = [I - pQ(*)] [I - P Q T (-A;)] (14) 

where Q T (— k) is the complex conjugate and transpose of Q(k). The first matrix is 
non-singular in the upper half complex A;-plane, while the second is non-singular in the 
lower half complex /c-plane. 

It can be shown that the factored correlation functions must be of the form 

Q(fc) = I - p / dre ikr Q(r) (15) 
where we used the following definition 

Aji = - °i) ( 16 ) 
S(r) = Q(r) - J dnQirJpQ^in - r) (17) 

Similarly, from Eq. (|T4*j) and Eq. (jHJ) we get, using the analytical properties of Q 
and Cauchy's theorem 

J(r) = Q(r) + y dr x J(r - rx)pQ(ri) (18) 

The general solution is discussed in [2312111211, and yields 

At 

mir) = <fi(r) + £ D$e-" X Jt < r (19) 

n=l 

g°(r) = (l/2)^[(r - ^/2) 2 - (^/2) 2 ] + /3,[(r - ^/2) - fa/2)] 

M 

+ E4 n)e " 2 " <Tj/2 [ e_2n(r_<Jj/2) - e_2 " ,Jl/2 ] A ii <r<ay (20) 

n=l 

The parameters Aj, (3j,Cij, are defined below and in Appendix 1. 
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3. The Laplace Transforms and the structure functions 

From Eq. (fT%|) we obtain the Laplace transform of the pair correlation function 

2vr E 9u(s) [S ej - ptqtj(is)] = (is) (21) 



where 



&(is)= / dre-^qUr)}' 



and 



where 



[(1 + sai/2)Aj + sPAe-^/s 2 - E -^—e^+^^C^f (22) 

m S + Z m 



' >A '/-,!'•<) = afM^i)Aj + a 2 Ms°i)f3j + (23) 

m S ~T Z m 



2 

Vi(ar) = [1 - x/2 - (1 + x/2) e -*]/(x 3 ) = — e^'H^zo jl) 
ii(x) = — (sinhx — x coshx) 



x 2 



[1 _ e -»| 

0i (x) = [1 - x - e^]/(x 2 ) = x^i(x) - o (x)/2; </> Q (x) = - 1 

x 

(24) 



In the factored case (see eq(jH)) 



4 m) = -4 m) 4 m) ; 4 m) = - ^- af»; = 2 7 rE^ir ) 4 m) 



(25) 



We remember that 



M 



m=l 



(26) 



and 



g y (is) = e sA ^of^i(s0-i)-4.j + e sAji of 0i(so-i)/3,- - 

5 | ^7~}~ 



m u 1 ^"i 



(27) 
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which can be written in the form 

My(s) = Stj - piqij(is) 



with 



M M+2 

8ij - dibj - Cidj -J2 e i m) fj m) = afh^ 



(2? 



a=l 



7T 

A 

2tt 



l0 - 2 Mso- t y^ 2 dj = p^ 2 = py s ^ 2 + 4 n)A(n) 



(m) / 



Pl e s ^ 2 a 2 V t\sa % ) 

e z m a t /2 f^M 



f 



O) = a ( m ) e ZmCTj/2 



„(m) 



l+^(l_ e -^ 

s 



Notice that eq. (j22|) can be written as 



with 



M M+2 
m a=l 



(l + SCTj/2) 

e 2 



e 2 



(30) 



~(m) _2i 

e- ; = e 2 



S + 2 r 



B 



(m) 



f (m) 



With this notation we can rewrite eq.([21|) as 

k 

After some algebra we get the general and simple result 



2Tvg ij (s) 



DJs) 



with 



L> = Det 



(31) 
(32) 



(33) 



(34) 
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We now apply the general result eq.(|33J) to a few selected examples of interest: 



8 



4-1. The Poly disperse hard spheres mixture 
For hard spheres 

<2 



J A 

Pi = 



2A 3 



2tt 



i + 1 ft 



Then we get, using the definitions of eq 

Ci = PiOifafai) dj = /3 C 



and also 



a, 

Ci 



(1 + sai/2)e- s ^/s 2 



The pair correlation function is 

2-KgiAs) -- 



s 2 D T 



(35) 



(36) 



4tt 2 



l> t = i -J2pi a j \ a jM s Vj) A °j - M s °j)Pj + -^tY. pi a fMs°i)M s °j)^ 



(37) 



4-2. The electrolyte limit 
We take the limit [27] 

of the general definitions. Then, 

Ci = piVifaisCTi) 



L n 0p - A . 



dj = $ + ajAjv = Pj 



p^o^Xso, 

i r 7T 



/ 



(0) 



afc-^l- 1 



(38) 



(39) 
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${is) = [(1 + aaifflAj + - £ 



S + Z r , 



-(s+z 



1 + so-t^Aj + s[3 3 } + sBfaf} 
{[(1 + 8^/2) A) + + s ^ (0) af } 



and also 



a, 



'1 + saj2) _«r t 
^ — -e 2 



e 2 
s 



s + z rt 



-(s)ai/2 e -z m ai/g(m) _ ^(m)> 



-8(Ti/2g(p) 



so that from eg. (|33jl we get the symmetric expression (compare |20| ) 



s 2 D, 



l + s 



with 



1 — {ab) —(ad) — (af) 
D± = Det -(cb) 1 - (cd) -(cf) 
-(e&) -(e^ 1 - (e/) 
where the scalar products ( see also j2D]) are defined by 

(ab) = <kbi, (cd) = £ dd h (ef) = s ^ j e i f i 

i i i 

and we use the definitions of eq. flBUjh 
5. THE PRIMARY CLOSURE 



The MSA closure condition obtained from Eq.(J7J) is 



i 

Using the results of the last section we got 



e 
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■jl a e 



E — rVEwW" } 



T (n) ^(m) 



= 
(46) 



6. The 1 Component Case: An explicit continued fraction solution 

In the one component case eg .([46)1 is simply 



with 



- 2nK n p [X n ] = z n f3 n 

5 n 



i+E 



m + ^rn 



-fa 



(47) 



X, 



(48) 
(49) 



L " j(n) + r n j( n ) 

fa = P&nX n 

In the diagonal approximation we get explicit equations for the scaling parameters (3. 
[T4] ( See appendix III., [II]) We recall that 

= a<t> {z n a) + pA J >^ 



and 



J (n) = 1 + pA J ' (n) 

A J -W = <7VoM^(7 - |W(^) 



(50) 
(51) 

(52) 
(53) 



and we have used the definition (j76j). 



The solution of the /? equations is obtained solving the linear equation 



(3 = M 



■2T 



(54) 



where 



M 



1 1 - 7i2 1 - 713 

1+712 1 1 - 723 

1 +713 1 + 723 1 





"A " 




' ri " 


? = 






r 2 
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2V 4- z — 2V — z 

'fnm ' 

We give the soultions for the first few cases [T""] : 



(i) 1 Yukawa (3 1 = 2T 1 

(ii) 2 Yukawas 

Pi 



02 



Ps 



2T 2 - 2I\ 



7l2 7l2 7l2 

(iii) 3 Yukawas: This case is slightly more complicated: 
The resolvent is 

( 723 ^ 
7? (3;1) = 731 , 

V 712 / 

so that the explicit solution is 

1 



Pi 

P2 = ~ 

Ps 



2T ■ 7? (3;1) 

723 I ; h S 2 3 I - ^23 



■ 1 ■ 






d 3 - 




713 ^ 



2T ■ ft 



(3;l) 



■ 1 " 






-d 3 - 




712 ^ 



d 3 

2t ■ 7i (3;1) 



+ Sl3 - ^13 



+ Si 2 - 212 



where d 3 = 712 + 723 - 7i3- 

2t ■ Jt^ 



Ps 



The general case of n > 3 can be found in the work of Blum et al. [Til 13*5] 
£.1. TTie 1 Yukawa case 

In the 1-Yukawa, 1 component case [27] the closure equation is simply 



2npK 1 [X 1 ] 2 = z 1 p 1 



1 + 



A 
2*i 



so that putting it all together we get the equation 



2/1 = 2r 1 (r 1 + z 1 ) [jw + r^ 1 



The physical branch yields the recursion relation 



Vi 



2npK 1 6 2 l 

Zl 



yi 



i(n+l) 



Vi 
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For large Z\ > 3, we may disregard the exponential terms in Eqs(|76|). and then we 
get the simple 



from where 



with 



ri 



2np (l + Y 



7T Z\ p 
4 



zfA< 



yi 



J 



(i) 



l 

Zl 



7rp(-2 + zi) 



z{ A 



2z 1 (l 



TT 2 p 2 _|_ 2 7T p (3 + 7T p) 
^2 A2 i 



2«f A 2 



3^ A 2 



2 i 



2(r 1 + z 1 )££ 1 



(64) 



(65) 



r.i 



1 + 



z[ A 



(ttV) 
2A 



+ 



27rp(l + ^) 2j/i (-A 



27rp(2-Zi) 

? 

-1 



2lA 



r) 7T 2 p 2 

Z ~~ IfA 2 " 



+ 



4?rp 1+ 



77 p 1 



z? A 2 



The solution is always convergent for the atractive case, and we do get the Onsage- 
rian limits for large densities and zero temperature correctly. It has been successfully 
tested numerically against other numerical methods [HZ] and Hernando, (unpublished, 
but using Pastore's criterion). For the repulsive case K > this equation has poles 
and unphysical divergencies occurring when the Onsagerian limits are attained. These 
divergencies are removed when a suitable reference hard core is introduced. 



6.2. The multiyukawa case 



The above solution is valid for any number of Yukawas n > 1. We write Eq. (|47|) in 
matrix form: 

1+AAn Pi/sis 
(3 2 /s 12 1+/Vs 2 2 • 



Vi 













01 

02 



(66) 



Using equations (|54*j) and (|B3j) we find 

1 + Pi/sn 0i/si2 

02/S12 1+/V S 22 



(67) 



with 



This relation is the extension of Eq. (j63|) to the multiyukawa case: The iteration follows 
exactly the same steps, only that now we have to use the extra equation (|54jl : The first 
iterate is 

1 

I 

1 



KO) 



J 212 1 £13 

S12 S13 

\ _|_ £12 ]_ \ _ £23 

512 S23 
\ _|_ £13 ^ I £23 

513 S23 



2/i/[^ (1) ] 2 



(68) 
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2 



(3 2 /si2 1+/V S 22 • 



(»>' 



where the superscript (n) indicates the level of the iteration. 



7. Summary of results 

• A new formulation of the pair distribution function for a large class of systems 
represented by a multiyukawa closure of the Ornstein-Zernike equation (Q) which 
includes systems such as water, ionic mixtures, ion channels and polymers is 
presented. This formulation is considerably simpler and more explicit than those 
in the literature because it involves directly measurable parameters such as the 
contact pair distribution function. It is also capable of handling the new improved 
versions of the MSA and the realistic octupolar model of water [7J. 

• The new result shown in equation (J2J) is simple and explicit: For example in the case 
of the SANS and SAXS experiments ^ |2] our result show explicitly the connection 
between the number of peaks in the diffraction spectra to the theoretical formula, 
because we get for dilute systems a sum of M Lorentzians, one for each yukawa 
term in the closure. 

• The one component case has been discussed in several papers, largely in 
collaboration with J. Hernando 6j, although most of the numerical results are 
still unpublished. The conclusion is that once the branch points are identified the 
continued fraction formalism always converges ( see also [37]). 



8. Appendix I: The algebraic solution of the general Yukawa closure of the 
Ornstein Zernike equation 



We quote the results from the review of Blum and Hernando [21]: The solution of the 
system of equations (|19|2()j) yields 



where 



3 A 



71 



1 + (1/2)^-^ 



(70) 

(71) 
(72) 

(73) 
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n (n) £(n) A (n) , 1 p(n). 

n i n) = 4 (n) + E^ n) ^ n) 

= 5jeaj(j)o(z n ai) - 2p e (3jafifji(z n a e ) 



where 



with 



g n) = -4EMl n) K + ^Q'K)] 



Here we define 



Cn = E Pk°k 
k 

A = 1 - ttC 3 /6 

/■OO 

9ij(s) = / drrgij(r)e 



lij = 27rg ij (z n )p j /z n 



1pl(s(Ji) = -[7 Z (3, SCTj) -7a (2, SCTj)] 



with the incomplete 7 function 



lz{n,z) 



{n-l)\ 



n-1 
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Mij-(s) = 6ij - Hij\ Hij = a i b . 



a=l 



(85) 



The determinant of this matrix can be expanded in its minors Aiy 

D T = Det\S i:j - Hij] = Y^[$ij ~ Vij] {-) i+j Mn ; Vj 

i 

This determinant can also be condensed into a simpler and more compact form 
which can also be expanded in its minors M^" 8 ) 

= ^ _ a («) . 6 (fl) ^(a,/?). .4(a,/J) = [_]«+/3 M (°./3) 

Using Cramer's rule then we find the general inverse 

[M-^j = 5 l3 + i- E ai a) 6f {^} ; V{a, 0} 



(86) 
(87) 



(89) 



J2{5 ik -ii ik }.[M-\s) 



\kj 



E 



(7)1,(7) 



— ^ij 



which means that we must satisfy 



7,0 



n 2^°k a k 



D 



(90) 
(91) 



E (Mm- 1 ^)]^ ) = E 4 Q) E *T [~ E 4 7) &f {A,*}) ; ^ 



1 E (_)^ a H 6 f }] + E a H (_ J_ E 



— A*ij 



(92) 



as it should 
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